clearvars -except Globaloption option
data = readtable('../../data/main_VAR/UK_source_pre1946.xlsx', 'Sheet', 'Full_Data_for_VAR');

date = data.Year;

startdatenum = 1729;
enddatenum = 1946;
middatenum = 1794;

startdate = find(date == startdatenum);
enddate = find(date == enddatenum);
middate = find(date == middatenum);
T_post = length(date(startdate:enddate));

trancate_year = 1850;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Load data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Tstart = startdate;
Tend = enddate;
T = Tend - Tstart + 1;

date = date(startdate:enddate);

data = data((data.Year >= startdatenum - 1), :);
data = data((data.Year <= enddatenum), :);

taxrevgdp = data.revtogdp(2:end); % government tax revenue to GDP ratio
taxrevgdp0 = data.revtogdp(1); % government tax revenue to GDP ratio
spendgdp = data.spendingtogdp(2:end); % government spending before interest exdp. to GDP ratio
spendgdp0 = data.spendingtogdp(1); % government spending before interest exp. to GDP ratio
surplusgdp = taxrevgdp - spendgdp; % primary surplus to gdp
rpcgdpgr = data.realGdpGrowth(2:end); % real gdp growth (in logs)
inflation = data.inflation(2:end); % growth in GDP price deflator (in logs)

ynom1 = data.short_r(2:end) / 100; % nominal yield on a 1-year Treasury bill, expressed per annum

ynom10 = data.rate_10year(2:end) / 100; % nominal yield on a 10-year Treasury bond, expressed per annum

dy = data.dy(2:end);
pdm = log(1 ./ dy); % log price/dividend ratio

%% convenience yield

cy_short1 = xlsread('../../data/supplement/UK_cony.xlsx', 'UK basis (ST)', 'ck5:ck45'); % 1873-1913
cy_short2 = xlsread('../../data/supplement/UK_cony.xlsx', 'UK basis (ST)', 'ck57:ck63'); % 1925-1931

cy_short = [mean(cy_short1) * ones(1873 - startdatenum + 1, 1); cy_short1; mean(cy_short2) * ones(1925 - 1914, 1); ...
                cy_short2; mean(cy_short2) * ones(enddatenum - 1931, 1)] / 100;

cy_long1 = xlsread('../../data/supplement/UK_cony.xlsx', 'UK basis (LT)', 'ck5:ck45'); % 1873-1913
cy_long2 = xlsread('../../data/supplement/UK_cony.xlsx', 'UK basis (LT)', 'ck57:ck63'); % 1925-1931

cy_long = [mean(cy_long1) * ones(1873 - startdatenum + 1, 1); cy_long1; mean(cy_long2) * ones(1925 - 1914, 1); ...
               cy_long2; mean(cy_long2) * ones(enddatenum - 1931, 1)] / 100;

% cy is zero before 1794;
cy_long = [zeros(middate - startdate + 2, 1); cy_long(3 + middate - startdate:end)];
cy = cy_long;

ynom1 = ynom1 + [zeros(middate - startdate + 1, 1); cy_short(3 + middate - startdate:end)];
ynom10 = ynom10 +cy_long(2:end);

ynom1 = log(ynom1 +1);
ynom10 = log(ynom10 +1);

divgrm = diff(log(data.dy) + data.stockindex);

gdebt = data.debttogdplevel;
% adjust for convenience yield

% calibrate cy/gdp ratio mean
cy_data = mean([(cy_long1) / 100 .* gdebt(1873 - startdatenum + 2:1913 - startdatenum + 2);
                (cy_long2) / 100 .* gdebt(find(date == 1925) + 1:find(date == 1931) + 1)]);
cy_data1 = mean((cy_long1) / 100 .* gdebt(1873 - startdatenum + 2:1913 - startdatenum + 2));
cy_data2 = mean((cy_long2) / 100 .* gdebt(find(date == 1925) + 1:find(date == 1931) + 1));

datew0 = (1728:1946)';
seig_rev = zeros(size(datew0));
seig_rev(1:find(datew0 == 1794)) = 0;
seig_rev(find(datew0 == 1795):find(datew0 == 1872)) = cy_data1;
seig_rev(find(datew0 == 1873):find(datew0 == 1913)) = gdebt(find(datew0 == 1873):find(datew0 == 1913)) .* (cy_long1 / 100);
seig_rev(find(datew0 == 1914):find(datew0 == 1924)) = cy_data2;
seig_rev(find(datew0 == 1925):find(datew0 == 1931)) = gdebt(find(datew0 == 1925):find(datew0 == 1931)) .* (cy_long2 / 100);
seig_rev(find(datew0 == 1932):find(datew0 == 1946)) = cy_data2;

taxrevgdp = taxrevgdp + seig_rev(2:end);
taxrevgdp0 = taxrevgdp0 + seig_rev(1); % use gdebt in year 1729 as the value in 1728;

deltalogg = [log(spendgdp(1)) - log(spendgdp0); log(spendgdp(2:end)) - log(spendgdp(1:end - 1))];
deltalogtau = [log(taxrevgdp(1)) - log(taxrevgdp0); log(taxrevgdp(2:end)) - log(taxrevgdp(1:end - 1))];

yieldspr = ynom10 - ynom1; % spread between 10-yr and 1-yr yields

%post 1850

taxrevgdp0 = taxrevgdp(trancate_year - startdatenum);
spendgdp0 = spendgdp(trancate_year - startdatenum);

rpcgdpgr = rpcgdpgr(trancate_year - startdatenum + 1:end);
ynom1 = ynom1(trancate_year - startdatenum + 1:end);
ynom10 = ynom10(trancate_year - startdatenum + 1:end);
yieldspr = yieldspr(trancate_year - startdatenum + 1:end);
inflation = inflation(trancate_year - startdatenum + 1:end);
deltalogtau = deltalogtau(trancate_year - startdatenum + 1:end);
deltalogg = deltalogg(trancate_year - startdatenum + 1:end);
taxrevgdp = taxrevgdp(trancate_year - startdatenum + 1:end);
spendgdp = spendgdp(trancate_year - startdatenum + 1:end);
surplusgdp = surplusgdp(trancate_year - startdatenum + 1:end);
divgrm = divgrm(trancate_year - startdatenum + 1:end);
pdm = pdm(trancate_year - startdatenum + 1:end);

T = enddatenum - trancate_year + 1;

y0nom_1 = mean(ynom1);
y0nom_10 = mean(ynom10);
yspr0 = mean(yieldspr);
pi0 = mean(inflation);
x0 = mean(rpcgdpgr);
tau0 = mean(deltalogtau);
g0 = mean(deltalogg);
surplus0 = mean(surplusgdp);

% log div/gdp ratio growth
divgrm = divgrm - inflation;

A0m = nanmean(pdm);
% mu_m = nanmean(divgrm);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Detrend tax and spending
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

tts(1) = taxrevgdp(1); 
gts(1) = spendgdp(1); 

for t = 2:T
    gts(t) = gts(t - 1) * exp(deltalogg(t) - g0);
    tts(t) = tts(t - 1) * exp(deltalogtau(t) - tau0);
end

%% dd colony government finance
% without India
data_colony_surplus = readtable('../../data/UK_colony/colony_surplus.xlsx');
date2 = data_colony_surplus.year;
samplestart = min(date2);
sampleend = max(date2);

rev_colony = data_colony_surplus.tau;
spending_colony = data_colony_surplus.g;
surplus_colony = rev_colony - spending_colony;

rev_colony = [zeros(samplestart - 1729, 1); rev_colony; zeros(1946 - sampleend, 1)];
spending_colony = [zeros(samplestart - 1729, 1); spending_colony; zeros(1946 - sampleend, 1)];

% India government finance data
Indiasurplusdata = readtable('../../data/UK_colony/India.xlsx');
samplestart = 1839;
sampleend = 1919;

India_Revenue = fillmissing(Indiasurplusdata.GROSSREVENUE_, 'constant', 0);
India_Spending = fillmissing(Indiasurplusdata.GROSSEXPENDITURE_, 'constant', 0);
India_Revenue = [zeros(samplestart - 1729, 1); India_Revenue; zeros(1945 - sampleend, 1)];
India_Spending = [zeros(samplestart - 1729, 1); India_Spending; zeros(1945 - sampleend, 1)];

year = [1729:1:1946]';

data_UK = xlsread('../../data/UK_colony/colonial_debt_cleaned.xls', 'UK');
INmat = xlsread('../../data/UK_colony/colonial_debt_cleaned.xls', 'IN');
INGDP = INmat(1:end, 3);
UKGDP = data_UK(:, 2);

India_Revenue_UKGDP = India_Revenue ./ UKGDP / 1000000;
India_Spending_UKGDP = India_Spending ./ UKGDP / 1000000;
India_Surplus_UKGDP = India_Revenue_UKGDP - India_Spending_UKGDP;

%combine
rev_colony = rev_colony(trancate_year - startdatenum + 1:end);
spending_colony = spending_colony(trancate_year - startdatenum + 1:end);
India_Revenue_UKGDP = India_Revenue_UKGDP(trancate_year - startdatenum + 1:end);
India_Spending_UKGDP = India_Spending_UKGDP(trancate_year - startdatenum + 1:end);

taxrevgdp = taxrevgdp + rev_colony + [India_Revenue_UKGDP; 0];
spendgdp = spendgdp + spending_colony + [India_Spending_UKGDP; 0];

%adjust for primary spending
UKmat = xlsread('../../data/UK_colony/colonial_debt_cleaned.xls', 'UK');
INmat = xlsread('../../data/UK_colony/colonial_debt_cleaned.xls', 'IN');
SAmat = xlsread('../../data/UK_colony/colonial_debt_cleaned.xls', 'SA');
AUmat = xlsread('../../data/UK_colony/colonial_debt_cleaned.xls', 'AU');
CAmat = xlsread('../../data/UK_colony/colonial_debt_cleaned.xls', 'CA');
NZmat = xlsread('../../data/UK_colony/colonial_debt_cleaned.xls', 'NZ');
colony_debt_gdp = (SAmat(:, 2) + AUmat(:, 2) +CAmat(:, 2) + INmat(:, 2) + NZmat(:, 2)) ./ UKmat(:, 2);
interest_expense_approximate = [colony_debt_gdp(trancate_year - startdatenum + 1:end); 0] .* ynom10;
spendgdp = spendgdp - interest_expense_approximate;

deltalogg = [log(spendgdp(1)) - log(spendgdp0); log(spendgdp(2:end)) - log(spendgdp(1:end - 1))];
deltalogtau = [log(taxrevgdp(1)) - log(taxrevgdp0); log(taxrevgdp(2:end)) - log(taxrevgdp(1:end - 1))];

y0nom_1 = mean(ynom1);
y0nom_10 = mean(ynom10);
yspr0 = mean(yieldspr);
pi0 = mean(inflation);
x0 = mean(rpcgdpgr);
tau0 = mean(deltalogtau);
g0 = mean(deltalogg);
surplus0 = mean(surplusgdp);

run ../../tools/run_var_estimate.m

save(['MAT/VARwithcolony_post1850.mat'],'-regexp', '^(?!(option|Globaloption)$).');

% CS computation 

%% Estimation
output_rp = 0.03 - mean(cy);
run ../../tools/cs_estimation;

%% steady-state upper bound calculation
upper = exp(pxbar) * (mean(taxrevgdp - spendgdp));

s_withcolony = s;
save ('MAT/colonyresult_post1850.mat', 's_withcolony');

